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1.  Six  examples  of  ionospheric  conditions  around  twilight.  Each  example  shows  the 

backseat  ter  power  (top  panel),  vertical  plasma  drift  (middle  panel),  and  zonal  plasma 
drift  (lower  panel),  plotted  against  the  scales  shown.  Symbols  plotted  in  the  middle 
panel  represent  the  vertical  drift  velocity  at  450 -km  altitude  for  reference.  The  left 
and  right  columns  portray  events  when  topside  spread  F  did  not  and  did  occur, 
respectively.  Note  that  UT  =  LT  +  5h .  3 

2.  Summary  Jicamarca  data  for  the  full  profile  analysis.  Prom  top  to  bottom,  the 
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03-06  LT  when  the  galaxy  was  directly  over  the  radar .  7 

3.  Example  occupation  profile  inversion  comparing  different  methods,  (a)  Simulated 
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1.  INTRODUCTION 


This  is  a  project  to  provide  mission  support  for  the  Communication/Navigation  Outage  Forecast 
System  (C/NOFS)  under  BAA  VS-03-01  during  its  first  four  years  of  operation.  The  objective  of  the 
overall  C/NOFS  satellite  mission  is  to  develop  and  test  methods  for  forecasting  so-called  equatorial 
spread  F  (ESF)  events  in  the  postsunset  ionosphere  associated  with  communications  interruptions  and  to 
collect  the  data  needed  to  do  so  on  an  operational  basis.  Cornell  was  to  support  the  mission  with  ground- 
based  radar  observations  of  background  ionospheric  parameters  and  of  ESF  events  from  the  Jicamarca 
Radio  Observatory  near  Lima,  Peru.  Jicamarca  is  capable  of  measuring  ionospheric  electric  field  and 
conductivity  profiles  from  the  valley  region  well  into  the  topside  ionosphere.  Jicamarca  can  also  make 
very  detailed  observations  of  ionospheric  irregularities  associated  with  ESF.  Three  main  tasks  were  to  be 
supported  under  this  project:  (1)  in-flight  calibration  of  the  electric  field/drift  meter  sensors  on  board  the 
spacecraft,  (2)  provision  of  Jicamarca  radar  support  for  experimental  campaigns,  both  before  the  launch 
and  during  the  satellite  operations,  and  (3)  analysis  and  interpretation  of  the  data  obtained  during  these 
campaigns  in  light  of  the  C/NOFS  mission  goals. 


2.  METHODS  AND  PROCEDURES 


The  C/NOFS  satellite  is  still  on  the  ground  awaiting  launch.  Consequently,  the  Cornell  project 
underwent  a  series  of  redefinitions  negotiated  with  our  science  counterparts  at  the  Air  Force  Research 
Laboratory  (AFRL).  Some  of  our  activities  involved  supporting  other  satellites,  including  DMSP  and 
COSMIC.  Many  of  the  original  objectives  therefore  carried  over  to  these  satellites.  Data  from  Jicamarca 
were  acquired,  processed,  and  analyzed  in  conjunction  with  numerous  DMSP  and  COSMIC  passes  during 
the  project  interval,  and  the  results  have  been  disseminated  to  the  C/NOFS  science  team  electronically 
and  through  meetings  and  workshops.  In  some  instances,  new  techniques  for  satellite  support  were 
developed.  Chief  among  these  is  a  mode  that  measures  calibrated  plasma  density  and  vector  electric  field 
profiles  simultaneously.  The  various  datasets  are  available  for  download  now  through  the  Madrigal  UAF 
server.  The  C/NOFS  satellite  will  receive  comparable  support  for  calibration,  validation,  and  science 
operations  upon  launch. 

Cornell  was  also  charged  with  delivering  ionospheric  measurements  to  initialize  and  validate 
climatological,  background,  and  forecast  models  under  development  at  AFRL.  These  data  were  taken, 
processed,  and  made  available.  Finally,  Cornell  was  to  investigate  ESF  variability  theoretically  and  to 
propose  new  forecast  strategies  based  on  our  findings.  This  component  of  the  project  also  proceeded  as 
planned. 

The  results  section  of  this  report  details  outcomes  from  three  specific  investigations.  These 
include  1)  developing  an  improved  forecast  strategy  for  ESF  based  on  our  evolving  understanding  of  a 
new  class  of  dynamical  plasma  instabilities  thought  to  initiate  the  phenomenon,  2)  improving  our  analysis 
of  Jicamarca  long-pulse  radar  data  to  better  support  the  radio  occultation  studies  being  performed  with 
COSMIC  as  well  as  general-purpose  DMSP,  C/NOFS,  and  model  calibration/validation  (CalVal),  and  3) 
investigating  satellite  radio  occultation  inversion  strategies  and  proposing  improvements  rooted  in 
contemporary  practices.  Each  of  these  outcomes  was  already  documented  in  separate  publications. 
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3.  FINDINGS 


Findings  pertaining  to  forecasting  ESF  based  on  a  new  instability  model  were  presented  in  [31]. 
A  new  method  of  analyzing  Jicamarca  long-pulse  radar  data  suitable  for  DMSP  and  C/NOFS 
CalVal  and  for  validating  COSMIC  occultation  measurements  was  described  by  [32].  Finally,  an 
improved  method  for  inverting  satellite  occultation  data  well  suited  for  E  and  lower  F  region  work 
was  promoted  by  [28].  The  most  important  findings  are  motivated  and  summarized  below.  Each 
of  the  following  sections  concludes  with  a  recommendation  for  the  ongoing  C/NOFS  mission. 


3.1.  Collisional  Shear  Instability  and  Equatorial  Spread  F 

A  striking  characteristic  of  the  bottomside  F  region  ionosphere  extending  from  late  afternoon  to  a 
few  hours  before  midnight  is  shear  flow.  Plasma  drifts  westward  in  the  bottomside  F  region  while 
drifting  eastward  near  and  above  the  F  peak,  with  the  shear  being  strongest  just  after  sunset.  That 
shear  flow  should  exist  in  the  bottomside  equatorial  F  region  was  evident  in  the  numerical  modeling 
results  of  [22],  argued  on  theoretical  grounds  by  [10],  and  observed  by  [37]  and  [59].  This  project 
focused  attention  on  shear  flow  in  the  postsunset  bottomside  F  region  ionosphere  and  the  possible 
role  it  plays  in  the  variability  of  equatorial  spread  F. 

While  the  existence  of  shear  flow  has  been  recognized  for  more  than  two  decades,  combined 
coherent  and  incoherent  scatter  observations  at  Jicamarca  allowed  us  to  locate  the  shear  node  and 
to  assess  the  shear  strength  quantitatively  for  the  first  time  (see  Fig.  1).  Shear  flow  exists  for 
a  period  of  several  hours  surrounding  sunset  and  is  most  intense  after  sunset  at  low  bottomside 
altitudes  approaching  the  valley  region.  The  scale  length  for  the  shear  is  measured  in  tens  of 
kilometers  there. 

Shear  flow  arises  mainly  from  the  competition  between  the  E  and  F  region  dynamos  and  from 
demands  on  the  vertical  boundary  currents  in  the  bottomside  imposed  from  above  by  the  imper¬ 
fectly  efficient  F  region  dynamo  and  from  below  by  interruptions  in  the  equatorial  electrojet  around 
the  evening  terminator  (20,  19].  Assessing  the  relative  importance  of  these  factors  quantitatively  is 
impeded  by  a  lack  of  neutral  wind  profile  measurements  in  the  equatorial  thermosphere  as  well  as 
of  off-equatorial  E  and  valley  region  conductivities. 

One  consequence  of  vertical  shear  is  the  existence  of  regions  where  the  zonal  plasma  flow  is 
strongly  retrograde,  i.e.,  in  the  direction  opposite  the  local  neutral  wind.  When  horizontal  con¬ 
ductivity  gradients  are  also  present  in  such  a  region,  wind-driven  interchange  instabilities  and 
plasma  irregularities  result.  This  now  appears  to  be  the  explanation  for  bottom-type  scattering 
layers,  which  occupy  retrograde  strata  and  appear  when  the  tilted  postsunset  F  region  supplies  the 
necessary  zonal  conductivity  gradient  [35]  (i.e.  see  the  narrow  layers  in  each  panel  of  Fig.  1) . 

On  occasions  when  topside  spread  F  occurs,  the  bottom-type  layers  seem  to  be  organized  in 
clusters  separated  by  tens  or  hundreds  of  km  [31].  We  know  this  from  interferometric/imaging 
experiments  at  Jicamarca  and  also  from  coherent  and  incoherent  scatter  radar  scans  from  the 
Altair  radar  on  Kwajalein  Atoll.  We  interpreted  this  modulation  as  evidence  of  intermediate 
and  large-scale  undulations  or  waves  in  the  bottomside  F  region.  The  wind-driven  instabilities 
would  be  stable  and  unstable  in  the  different  phases  of  such  waves.  In  most  instances,  the  topside 
irregularities  that  followed  also  exhibited  structuring  with  the  same  characteristic  scale  size.  The 
intermediate  and  large-scale  waves  therefore  appear  to  serve  as  seed  waves  for  full-blown  spread  F. 
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Figure  1:  Six  examples  of  ionospheric  conditions  around  twilight.  Each  example  shows  the  backscat- 
ter  power  (top  panel),  vertical  plasma  drift  (middle  panel),  and  zonal  plasma  drift  (lower  panel), 
plotted  against  the  scales  shown.  Symbols  plotted  in  the  middle  panel  represent  the  vertical  drift 
velocity  at  450  -km  altitude  for  reference.  The  left  and  right  columns  portray  events  when  topside 
spread  F  did  not  and  did  occur,  respectively.  Note  that  UT  =  LT  +  5h. 
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The  existence  of  seed  waves  is  often  postulated  to  explain  the  rapid  appearance  and  development 
of  radar  plumes  and  intense  plasma  irregularities  shortly  after  sunset  in  the  equatorial  ionosphere. 
Gravity  waves  propagating  into  the  thermosphere  may  seed  the  ionosphere  as  well  but  axe  not 
expected  to  preferentially  excite  waves  at  the  observed  wavelengths. 

Beyond  creating  conditions  for  wind-driven  instabilities,  shear  flow  may  itself  destabilize  the 
bottomside  F  region.  Boundary  value  analyses  by  [18],  [26],  and  [53]  indicated  that  shear  flow 
generally  stabilizes  ionospheric  interchange  instabilities,  which  was  the  prevailing  sentiment  for 
many  years.  However,  [16],  [50],  and  [15]  challenged  this  idea  by  pointing  out  the  limitations 
of  boundary  value  analysis,  which  is  incomplete  when  applied  to  non-normal  models  like  those 
describing  sheared  flows. 

By  extending  the  work  of  [34]  to  the  equatorial  zone,  [30]  showed  that  a  collisional  shear  insta¬ 
bility  may  exist  in  the  postsunset  ionosphere  and  be  responsible  for  the  precursor  waves  we  observe. 
In  the  collisional  limit,  retrograde  plasma  motion  rather  than  shear  per  se  drives  the  instability. 
However,  it  is  not  just  an  obvious  variant  of  the  standard  wind-driven  interchange  instability,  since 
instability  occurs  when  the  neutral  wind  and  the  background  density  gradient  are  normal  to  one 
another.  The  instability  is  the  collisional  branch  of  electrostatic  Kelvin  Helmholtz  and  is  therefore 
aptly  termed  a  collisional  shear  instability.  The  instability  is  not  damped  by  collisions  and  does 
not  require  a  prereversal  enhancement  and  the  elevation  of  the  Flayer  to  function.  This  may  allow 
it  to  operate  prior  to  the  onset  of  gravitational  interchange  instabilities,  which  are  stabilized  by 
collisions,  and  thereby  to  serve  as  a  source  of  seed  waves. 

An  approximate  linear  growth  rate  expression  for  the  collisional  shear  instability  was  derived 
in  magnetic  coordinates  (p,<?,0)  by  [31]: 


~  «o)n0') 

“n°>  +  ^  ^ 

where  the  angle  brackets  represent  averages  over  a  magnetic  flux  tube  as  well  as  over  a  range  of  apex 
altitudes,  the  prime  represents  a  derivative  with  respect  to  the  p  coordinate,  and  where  the  k.  terms 
are  constant,  nondimensional  wavenumbers  related  to  the  wavenumber  through  the  appropriate 
geometric  scale  factor,  e.g.  The  inherently  nonlocal  nature  of  the  instability  prevents 

the  derivation  of  an  exact  expression,  but  the  approximation  reveals  at  least  how  the  growth 
rate  scales  with  background  parameters.  Remarkably,  the  possible  range  of  growth  rated  greatly 
exceeds  that  for  the  generalized  Rayleigh  Taylor  instability,  highlighting  the  importance  of  this 
mechanism  in  the  initiation  of  ESF.  Retrograde  plasma  motion  and  a  steep  vertical  plasma  density 
gradient  in  the  retrograde  stratum  are  conducive  to  wave  growth.  Large  flux-tube  integrated  plasma 
conductivity  suppresses  growth.  Whether  the  former  overcomes  the  latter  depends  on  the  strength 
of  the  neutral  wind,  the  cause  of  retrograde  drifts,  and  the  off-equatorial  conductivity  distribution. 

For  example,  if  the  boundary  current  condition  is  mainly  responsible  for  retrograde  drifts  in 
the  bottomside,  then  the  bottomside  flux  tubes  need  not  be  heavily  loaded  by  conductivity  in  the 
off-equatorial  E  and  valley  regions,  and  shear  instability  could  be  robust.  If  shear  flow  is  mainly 
induced  by  a  strong  E  region  dynamo  controlling  flux  tubes  in  the  bottomside,  however,  then  the 
shear  instability  would  likely  be  damped.  The  existence  of  irregularities  in  the  bottom-type  layers 
gives  no  clue  about  which  scenario  is  operating,  since  they  can  exist  in  either  case. 

Thermospheric  wind  measurements  made  by  the  C/NOFS  satellite  will  be  critical  for  quantifying  the  shear 
instability  theory.  In  order  to  evaluate  (1),  information  about  the  off-equatorial  conductivity  distribution  will  also  be 
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required.  This  could  be  provided  by  a  meridional  chain  of  sounders,  a  bistatic  HF  sounder  link,  or  by 
satellite  occultation  experiments.  The  same  information  is  required  both  to  quantify  the  causes  of  shear 
flow  in  the  equatorial  ionosphere  and  to  evaluate  the  growth  rate  of  the  gravitational  interchange 
instabilities  responsible  for  spread  F. 

3.2.  Jicamarca  Incoherent  Scatter  Data  Analysis  Using  a  New  Full-Profile  Approach 

The  protracted  pre-launch  period  for  C/NOFS  has  afforded  the  opportunity  to  reexamine  the 
analysis  of  Jicamarca  incoherent  scatter  data  used  to  support  spacecraft  missions.  Incoherent  scatter  data 
from  a  hybrid  long-pulse/double-pulse  experiment  run  at  Jicamarca  during  C/NOFS  and  COSMIC 
preparatory  campaigns  have  been  analyzed  using  a  full-profile  analysis  stemming  from  the  one  described 
by  [24],  In  this  case,  plasma  density,  electron  and  ion  temperatures,  and  light  ion  composition  profiles  in 
the  topside  were  estimated  simultaneously.  Full-profile  analysis  will  be  crucial  for  C/NOFS  support  at 
Jicamarca,  since  the  long  correlation  time  of  the  incoherent  scatter  signal  at  50  MHz  invalidates 
conventional  gated  analysis  of  the  long-pulse  data  we  plan  to  collect. 

For  most  of  its  history,  Jicamarca  employed  a  double  pulse  technique  for  incoherent  scatter 
experiments  [8,  48],  Pairs  of  short  pulses  with  opposite  circular  polarizations  were  transmitted,  and  lag 
product  matrices  were  assembled  from  the  backscatter,  one  lag  at  a  time.  The  size  of  the  scattering  volume  was 
set  by  the  pulse  lengths  rather  than  the  lag  spacing,  and  problems  associated  with  inhomogeneity  were  thereby 
avoided.  The  double  pulse  experiment  was  also  relatively  immune  to  clutter  from  the  ground,  the  equatorial 
electrojet,  and  other  sources.  Faraday  rotation  measurements  using  phase  information  from  the  orthogonal 
polarizations  yielded  an  independent  measurement  of  absolute  electron  density  for  calibration  purposes  [7], 
However,  the  double  pulse  experiment  was  very  slow  and  did  not  frilly  exploit  the  average  power  limitations 
of  the  transmitters.  Given  reasonable  incoherent  integration  times,  the  experiment  could  not  yield  useful 
ionospheric  profiles  at  protonospheric  heights. 


More  recently,  a  hybrid  experiment  combining  the  double  pulses  with  an  alternating  code  was 
implemented  at  Jicamarca  [27],  The  alternating  code  permits  the  estimation  of  the  lag  product  matrix  with 
range  resolution  set  by  the  bit  length  of  the  code  rather  than  the  pulse  length  [41,  42,  43],  The  range  resolution 
of  all  the  lags  is  uniform  except  for  the  zero  lag.  The  alternating  code  thus  resolves  many  analysis  problems 
associated  with  inhomogeneity.  The  alternating  code  is  also  much  faster  than  the  double  pulse,  providing 
multiple  lag  product  estimates  from  each  pulse  and  fully  utilizing  the  duty  cycles  of  the  transmitters.  Given 
reasonable  incoherent  integration  times,  the  hybrid  experiment  yields  useful  ionospheric  profiles  well  into  the 
protonosphere. 

While  representing  an  improvement  over  the  double  pulse  at  altitudes  where  it  can  be  applied,  the 
alternating  code  is  still  less  sensitive  than  the  long  pulse.  This  is  because  most  of  the  signal  in  a  given  receiver 
sample  comes  from  ranges  other  than  the  range  under  analysis.  The  experiment  is  self  clutter  prone,  with  a 
signal-to-clutter  ratio  never  exceeding  the  reciprocal  of  the  code  length.  Added  clutter  is  the  price  paid  for 
reduced  range  ambiguity. 

We,  therefore,  pursued  a  new  strategy  for  long-pulse  analysis  introduced  by  [24]  and  investigated 
by  [44].  Rather  than  manipulating  the  lag  products  from  a  long-pulse  experiment,  constructing 
approximations  of  autocorrelation  functions,  and  fitting  the  results  for  ionospheric  parameters  one  range 
at  a  time  (so-called  “gated”  analysis),  [24]  proposed  the  simultaneous  estimation  of  complete  parameter 
profiles  directly  from  the  lag  product  matrices  as  they  are  measured. 
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The  technique  is  termed  “full  profile  analysis”  and  belongs  to  a  class  of  problems  in  statistical 
inverse  theory.  The  technique  neither  degrades  the  data  using  summation  rules  nor  introduces  clutter  with 
pulse  coding.  In  practice,  it  entails  solving  the  forward  model  iteratively  until  a  cost  function  composed 
of  the  model  prediction  error  (essentially,  the  chi-squared  statistic)  and  a  combination  of  model  metrics  is 
minimized. 

The  theory  that  relates  the  autocorrelation  function  p(k,r;r)  to  the  state  variables  of  the 
ionospheric  plasma  (electron  density,  electron  and  ion  temperatures,  ion  composition,  drifts,  etc.) 
is  well  known  and  need  not  be  reviewed  here  (see  for  example  [4,  9,  5,  57,  11,  12,  13,  51,  52,  23,  60]). 
Application  of  the  theory  at  Jicamarca  is  complicated  by  the  small  magnetic  dip  angle,  however.  A 
considerable  body  of  work  exploring  the  theory  at  small  aspect  angles  has  emerged  in  recent  years. 
Here,  we  highlight  findings  relevant  to  the  forward  model. 

[36]  first  noted  that  the  incoherent  scatter  spectrum  at  small  aspect  angles  is  even  smaller  than 
the  standard  theory  predicts.  (See  also  [38].)  At  the  same  time,  it  had  long  been  known  that 
the  standard  theory  failed  to  predict  the  correct  electron-ion  temperature  ratio  for  finite  but  small 
aspect  angles,  of  the  order  of  a  few  degrees  or  less  [48].  [56]  were  able  to  account  for  both  of  these 
phenomena  (qualitatively  for  the  former  and  quantitatively  for  the  latter)  by  including  the  effects 
of  electron  Coulomb  collisions  in  the  calculation  of  the  electron  admittance.  Their  calculations 
were  based  on  Monte  Carlo  simulations.  Summarizing  these  findings,  [1]  went  on  to  show  that 
the  new  Coulomb  collision  corrections  bring  temperature  ratio  observations  made  a  few  degrees 
off  perpendicular  at  Jicamarca  into  compliance  with  expectation.  Those  corrections,  which  have 
been  made  available  in  tabular  form  (M.  Sulzer,  personal  communication)  axe  incorporated  in  the 
present  analysis. 

A  second  issue  involves  the  relationship  between  plasma  density,  the  radar  scattering  cross- 
section,  and  power  profile  measurements.  At  small  aspect  angles,  the  incoherent  scatter  cross 
section  depends  on  the  electron  and  ion  temperatures  and  the  plasma  density  in  a  complicated 
way  not  predicted  by  the  approximation  for  large  aspect  angles  given  by  [6]  employed  widely  at 
other  ISR  facilities.  [46]  have  examined  the  effect  in  detail,  including  for  the  first,  time  the  effects 
of  electron  Coulomb  collisions  in  their  calculations.  Extending  the  work  by  [61],  they  incorporate 
the  effects  of  Coulomb  collisions  analytically  with  the  introduction  of  a  heuristic  operator  in  the 
so-called  Gordeyev  integral  which,  when  evaluated  numerically,  reproduces  the  results  of  [56].  [49] 
recently  used  their  results  to  convert  Jicamarca  power  profiles  to  density  profiles,  finding  good 
agreement  with  Faraday  rotation  measurements.  We  adopt  the  same  methodology  here,  making 
use  of  Faraday  rotation  only  for  absolute  density  calibration. 

Full  profile  analysis  is  an  example  of  a  discrete  inverse  problem  of  the  form  G(m)  —  d,  where 
d  is  a  column  vector  representing  known  data,  m  is  a  column  vector  representing  unknown  model 
parameters,  and  G  is  the  known  theory  relating  them.  In  this  instance,  the  data  axe  lag  profile 
estimates,  the  model  parameters  are  plasma  density,  temperature  and  composition  specified  at 
various  altitudes,  and  the  theory  is  incoherent  scatter  theory  combined  with  the  radar  ambiguity 
function.  If  the  model  is  linear  or  its  estimate  is  sufficiently  close  to  the  solution  that  the 
theory  can  be  linearized  about  mest,  then  the  problem  becomes  simply  Gmest  =  d,  where  d  G  Rn, 
m  G  Rm,  and  G  G  Rnxm.  The  apparent  simplicity  of  even  this  problem  is  deceptive,  and  its 
solution  in  the  presence  of  noisy  data  almost  inevitably  necessitates  recourse  to  statistical  inverse 
methods.  For  an  introduction  to  inverse  methodology,  consult  for  example  [45],  [58],  and  [2]. 
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Figure  2:  Summary  Jicamarca  data  for  the  full  profile  analysis.  From  top  to  bottom,  the  panels 
represent  electron  density,  electron  temperature,  ion  temperature,  hydrogen  ion  fraction,  and  helium 
ion  fraction  (25%  full  scale).  A  data  gap  exists  between  03-06  LT  when  the  galaxy  was  directly 
over  the  radar. 
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Zero  lag  inversion  In  the  present  problem,  each  of  the  measured  long-pulse  lag  products  rep¬ 
resents  the  two  dimensional  convolution  of  the  radar  range-time  ambiguity  function  and  the  auto¬ 
correlation  function  of  the  scattering  medium.  The  zero  lag  is  no  exception.  However,  since  the 
ACF  is  nearly  uniform  over  the  small  delays  r  entering  into  the  problem,  we  can  regard  the  zero 
lag  as  the  one-dimensional  correlation  of  a  ‘true’  power  profile  with  the  emitted  radar  pulse  shape 
to  a  good  approximation.  Here,  ‘true’  means  the  profile  that  would  be  measured  ideally  if  short 
pulses  were  used.  We  can  therefore  recover  the  former  by  deconvolving  the  latter  from  the  zero  lag 
profile. 

Deconvolution  is  performed  using  Tikhonov  regularization.  This  is  accomplished  with  aug¬ 
mented  linear  least  squares  minimization,  the  fourth  enumerated  strategy  above.  Minimization 
is  performed  using  a  non- negative  least  squares  (NNLS)  routine  that  prevents  the  occurrence  of 
negative  results  anywhere  in  the  deconvolved  profile  [40].  This  strategy  enhances  the  stability  and 
robustness  of  the  methods  and  permits  smaller  regularization  parameters  than  would  otherwise 
be  possible.  It  is  noteworthy  that  the  droop  in  transmitter  power  that  occurs  during  long  pulse 
transmission  produces  a  forward  model  that  is  better  conditioned  than  would  be  a  uniform  pulse. 

The  ‘true’  power  profile  combined  with  incoherent  scatter  theory  and  the  shape  of  the  radar  pulse 
complete  the  forward  model.  Given  specifications  for  the  ionospheric  density,  temperature,  and 
composition  profiles  along  with  the  radar  range- time  ambiguity  function,  we  can  predict  the  absolute 
lag  product  matrix.  In  practice,  the  integrals  in  the  ambiguity  function  are  evaluated  discretely 
using  a  trapezoidal  rule.  The  proxy  electron  density  profile  is  a  sufficiently  good  approximation  for 
ne  forward  modeling  intents  and  purposes  and  need  only  be  corrected  at  the  conclusion  of 

the  analysis. 


Inverting  the  entire  lag  profile  matrix  For  the  forward  model,  electron  temperature,  the 
electron/ ion  temperature  ratio,  and  the  hydrogen  and  helium  ion  fractions  are  specified  coarsely  at 
20  uniformly-spaced  ranges  spanning  the  solution  space.  These  become  the  knots  used  for  cubic  B- 
spline  interpolation  [3].  The  interpolated  functions  thus  created  are  twice-differentiable  continuous 
and  provide  a  means  of  specifying  the  ionospheric  parameter  profiles  with  arbitrarily  fine  resolution 
for  forward  modeling.  Computing  the  B-spline  coefficients  is  fast  and  efficient,  and  experimentation 
has  shown  that  the  topside  profiles  measured  with  conventional  gated  analysis  at  Jicamarca  can  be 
well  represented  with  20  or  fewer  knots. 

The  model  parameters,  the  4x20  knot  values,  are  calculated  using  augmented  nonlinear  least 
squares  minimization.  The  minimization  algorithm  is  conventional  Levenberg-Marquardt.  The 
parameters  are  initialized  with  values  resulting  from  an  extensive  grid  search  when  a  large  family 
of  representative  trial  profiles  are  generated  and  ax2  statistic  calculated  for  each.  It  is  generally 
possible  to  find  initial  profiles  this  way  for  which  x2 /m  <  5,  m  being  the  number  of  parameters. 
Levenberg  Marquardt  iteration  generally  yields  final  values  for  x2/m  ~  2-3.  Smaller  values  can  be 
produced  only  at  the  cost  of  producing  very  irregular  parameter  profiles. 

Augmentation  in  this  case  means  introducing  six  additional  cost  functions  that  help  to  reg¬ 
ularize  the  output  of  the  algorithm.  Each  of  these  is  multiplied  by  an  adjustable  weight.  The 
weights  have  been  set  so  as  to  inhibit  results  that  are  inconsistent  with  general  expectations  for  the 
ionosphere.  The  cost  functions  are  enumerated  below. 

1.  (\d2Te/dz2\) 
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2.  (\d?Ti/dz2\) 

3.  {\d2H+/dz2\) 

4.  exp(—  min(0, 1  —  H+  —  He+)) 

5.  (He+) 

6.  exp(-mm(0,Te-Ti)) 

The  first  three  coat  functions  pertain  to  the  global  range  averages  of  the  second  derivatives  of 
the  electron  temperature,  ion  temperature,  and  hydrogen  ion  fractions,  respectively.  Suppressing 
these  quantities  inhibits  spurious  oscillations  in  the  profiles  associated  with  poor  conditioning.  Note 
that  the  net  hydrogen  fraction  in  our  algorithm  is  defined  by 

H+  =  (1  +  tanh  ff+)/2 

A  similar  function  is  used  to  express  the  helium  ion  fraction.  The  functions  restrict  the  ion  fractions 
to  values  between  0  and  1.  It  remains  possible  for  the  combined  hydrogen  and  helium  ion  fractions 
to  exceed  unity,  and  so  a  fourth  cost  function  was  introduced  to  inhibit  this. 

A  fifth  cost  function  is  a  very  weak  prohibition  of  helium  ions  overall.  This  is  meant  to  prevent 
the  algorithm  from  producing  spurious  He+  layers  that  might  otherwise  come  about  as  a  conse¬ 
quence  of  any  systemic  errors  in  data  acquisition  and  processing.  Helium  ions  indicated  by  the  full 
profile  analysis  must  be  pointed  to  by  clear  signatures  in  the  lag  product  data.  Finally,  a  sixth 
cost  function  is  a  prohibition  against  electron-to-ion  temperature  ratios  less  than  unity.  This  is  a 
condition  that  we  expect  never  to  occur  in  nature,  and  the  cost  function  prevents  the  full  profile 
algorithm  from  exploring  that  region  of  solution  space. 

Twenty-four  hours  of  hybrid  double/long-pulse  data  were  taken  during  April  19-20,  2006,  at 
Jicamarca  and  processed  using  the  full-profile  algorithm.  In  the  analysis,  only  the  diagonal  compo¬ 
nents  of  the  error  covariance  matrices  were  used.  While  the  formalism  surrounding  the  calculation 
and  incorporation  of  the  complete  covariance  matrix  is  straightforward,  the  storage  and  computa¬ 
tional  demands  are  significant,  and  sparse  matrix  solvers  are  required  in  practice.  Sparse  matrix 
math  has  yet  to  be  incorporated  in  our  algorithm. 

Processed  data  appear  in  the  summary  plot  shown  in  Fig.  2  for  the  24-hour  run.  Data  between  0300- 
0600  LT  when  the  galaxy  was  inside  the  antenna  beam  have  been  excised,  as  have  intervals  when 
residual  echoes  from  space  debris  interfered  with  the  analysis.  While  there  is  significant  scatter  in  some  of 
the  parameters,  particularly  in  the  postsunrise  temperatures,  this  is  a  consequence  of  low  signal  levels  and 
consistent  with  our  error  analysis. 

The  salient  features  of  Fig.  2  are  the  dramatic  electron  heating  after  sunrise,  the  low  and  rela¬ 
tively  constant  0+ /H+  transition  height,  and  the  appearance  of  a  thin,  patchy  helium  layer  in  the 
morning.  They  are  remarkably  similar  to  model  predictions  for  equinox  solar  minimum  conditions 
over  Jicamarca  made  by  [25]  except  that  the  measured  plasma  temperatures  are  somewhat  smaller 
at  sunrise.  These  features  have  also  been  suggested  by  gated  analyses  in  the  past  and  so  do  not 
appear  to  be  creations  of  the  full-profile  analysis. 

Full-profile  analysis  is  now  available  for  C/NOFS  CalVal.  Furthermore,  it  can  and  should  be 
used  as  a  baseline  against  which  to  tune  and  validate  model  ionospheres  such  as  the  one  under 
development  at  AFRL.  The  full-profile  results  (density,  temperature,  and  composition  profiles)  are 
the  most  reliable  representation  of  the  background  equatorial  ionosphere  available  to  date. 
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3.3.  Optimal  Radio  Occultation  Profile  Inversion 


Radio  occultation  measurements  of  the  Earth’s  atmosphere  and  ionosphere  using  global  positioning 
satellites  and  receivers  placed  on  satellites  in  low  Earth  orbit  (LEO)  are  an  important  new  source 
of  remote  sensing  information,  supplementing  ground-based  and  in  situ  measurements  and  serving 
as  sources  for  assimilative  models.  The  COSMIC  mission  in  particular  is  dedicated  to  specifying 
the  ionosphere  globally  for  purposes  of  assimilative  modeling  and  stability  forecasting.  Developed 
for  planetary  atmospheric  research  (e.g.,  [47, 14]),  the  radio  occultation  technique  involves  inferring 
gradients  in  the  index  of  refraction  from  differential  phase  measurements  and  the  implied  optical 
path  length  between  transmitting  and  receiving  satellites.  The  index  of  refraction  in  the  Earth  s 
atmosphere  is  a  function  of  temperature,  pressure,  water  and  water  vapor  content,  and  electron 
density,  where  the  last  factor  has  by  far  the  most  significant  effect.  Inferring  total  electron  content 
(TEC)  on  the  ray  path  between  the  satellites  is  then  a  relatively  straightforward  task.  Reviews  of 
different  analysis  techniques  have  been  given  by  [39]  and  [54]. 

The  purpose  of  this  project  was  to  investigate  a  pathology  that  arises  in  the  conventional  inver¬ 
sion  strategy  and  explore  means  of  mitigation.  Two  conceptually  related  methodologies,  Tikhonov 
regularization  and  maximum  entropy,  were  described  and  compared.  While  both  offer  substantial 
improvement  over  straight  linear  inverse  transformation,  the  latter  was  seen  to  outperform  the 
former  in  its  ability  to  retain  fine  structure  while  suppressing  artifacts. 

We  are  interested  in  the  specific  problem  of  recovering  altitude  profiles  of  ionospheric  electric 
density  from  radio  occultation  observations  made  with  satellite  pairs.  We  assume  spherical  sym¬ 
metry  so  that  the  electron  density  is  a  function  of  altitude  r  alone.  We  further  neglect  ray  bending 
and  birefringence  in  this  analysis  and  consider  only  the  effect  that  the  ionospheric  plasma  has  on 
straight-line  propagation  through  variations  in  the  index  of  refraction  n : 

n  =  1-80.48^ 

r 

where  ne  is  the  electron  number  density  in  MKS  units  and  /  is  the  radio  frequency  in  Hertz,  which 
is  taken  to  be  well  above  the  critical  frequency  in  the  ionosphere.  Because  of  the  dispersion,  signals 
received  at  different  frequencies  will  have  different  optical  path  lengths  and  arrive  out  of  phase.  For 
a  given  propagation  path  with  closest  approach  (to  the  center  of  the  Earth)  s,  the  phase  difference 
will  be  related  to  the  electron  density  profile  through  Abel’s  transform: 


*(«) 


ne(r) 


rdr 

yV2  —  s2 


(2) 


where  it  is  assumed  that  both  the  transmission  and  reception  points  are  above  the  altitudes  where 
significant  electron  density  exists.  (This  will  not  generally  be  true  for  LEO  satellites,  but  a  modified 
methodology  exists  in  that  case  that  uses  positive  and  negative  elevation  rays  to  resolve  the  problem 
[21].)  Here,  the  constant  C  for  the  two  frequencies  f\  and  ft  is  27t(/i  -  h)/{cf  1/2)-  (i.e.  not  the 
speed  of  light).  As  the  satellites  move,  they  sweep  out  the  continuous  curves  4>(s).  The  task  is  to 
recover  the  electron  density  profile  ne(r)  by  inverting  (2). 

In  the  event  that  ne(r)  decreases  more  rapidly  than  r-1  at  the  upper  bound,  the  inverse  of  (2) 
exists: 


n  (r) _ l  r  dA  dJ— 

Mr)  ~  nc  Jr  ds^rzr^ 


(3) 
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This  transformation  has  the  added  benefit  of  acting  on  the  derivative  of  the  phase  term,  sidestepping 
ambiguities  associated  with  absolute  phase  measurement.  It  is  a  straightforward  matter  to  evaluate 
(3)  using  standard  numerical  integration  techniques.  We  will  see  that  the  problem  is  well  posed 
and  that  (3)  implies  no  systematic  distortion  of  the  desired  profile.  However,  (3)  possesses  a  certain 
pathology  that  makes  it  poorly  suited  for  ionospheric  investigation.  That  pathology  and  a  remedy 
are  described  below. 

Consider  the  discrete  versions  of  (2)  and  (3).  The  electron  density  and  phase  difference  profiles 
can  be  discretized  spatially  and  regarded  as  column  vectors.  We  designate  the  phases  with  the 
vector  d  for  data  and  the  electron  densities  with  the  vector  m  for  model.  Let  d  be  specified  at  n 
points  and  m  at  m  points,  where  the  integers  n  and  m  need  not  be  equal.  The  integrals  can  then 
be  carried  out  using  summation  rules.  A  simple  trapezoidal  integration  scheme  will  be  used  here, 
although  higher  order  schemes  are  certainly  possible.  The  derivative  in  (3)  can  be  evaluated  with 
backward  differencing.  Finally,  the  forward  and  inverse  transforms  assume  the  form 

d  =  Gm  (4) 

m  =  G*d  (5) 

where  G  is  the  nxm  linear  transformation  implied  by  the  integral  in  (2)  and  G*  is  its  mxn  general¬ 
ized  inverse.  The  caret  in  (5)  denotes  an  estimate  with  accuracy  depending  on  the  qualities  of  G*; 
note  that  G*  need  not  and  will  not  in  general  be  G"1.  The  way  forward  is  the  subject  of  inverse 
theory,  which  is  described  in  a  number  of  standard  texts  (see  for  example  [2],  [45],  and  [58]).  Using 
the  Moore-Penrose  generalized  inverse  and  given 

G  =  UpSpVj,  (6) 

&  =  vps;%  (7) 

where  C/p,  Sp,  and  Vp  are  matrices  found  through  the  singular  value  decomposition  of  G,  truncated 
to  dimensions  limited  by  the  rank  p  of  G.  Here,  Up  and  Vp  are  orthogonal  matrices  (size  nxp  and 
mxp,  respectively)  spanning  the  data  and  model  ranges,  respectively,  and  Sp  is  the  pxp  diagonal 
matrix  composed  of  positive  definite  singular  values. 

For  the  purposes  of  this  analysis,  it  will  suffice  to  consider  the  even  determined  case  in  which 
G  is  full  rank  and  p  =  m  =  n.  Then  the  singular  values  are  all  nonzero,  VpVp  =  UpUp  =  /,  and 

Gf  =  VS-'U* 

=  {USV1)'1 

=  G'1  (8) 

The  inverse  is  unique  and  exact,  since  m  =  G^Gm  =  m. 

The  accuracy  of  the  model  estimate  in  the  presence  of  noise  depends  on  the  stability  of  the 
inversion  and  is  the  issue  motivating  this  analysis.  The  structure  of  (8)  gives  an  indication.  Here,  the 
inverse  operator  can  be  seen  to  decompose  the  data  vector  into  its  component  singular  vectors  and 
to  scale  each  one  by  the  corresponding  inverse  singular  value.  Ill-conditioned  problems  are  signaled 
by  ratios  of  the  smallest-to-largest  singular  values  competing  with  the  floating  point  precision  of 
the  computer,  since  roundoff  error  prohibits  accurate  retention  of  all  the  singular  modes,  leading 
to  instability.  For  the  Abel  inversion,  this  ratio  is  just  comparable  to  n.  Even  so,  this  suggests  that 
error  amplification  will  be  substantially  different  for  different  kinds  of  features  in  the  data.  This 
phenomenon  can  be  explored  further  through  an  example. 


II 


Fig.  3  illustrates  the  performance  of  the  Abel  transform  inversion  by  applying  the  linear  trans¬ 
formation  in  (5)  to  simulated  data.  For  these  simulations,  n  =  m  =  60.  Panel  (a)  shows  the  initial 
electron  density  profile,  n€(r),  discretized  on  a  grid  with  unit  spacing.  The  electron  density  here 
is  dimensionless  and  normalized  such  that  the  area  under  the  curve  is  also  unity.  The  profile  is 
composed  of  three  Maxwellian  components  and  is  intended  to  vaguely  resemble  the  ionosphere. 

Panel  (b)  shows  simulated  occultation  data  0(s)  computed  using  the  transformation  in  (4).  For 
these  calculations,  we  take  the  constant  C  in  (2)  and  (3)  to  be  unity  for  simplicity.  Two  curves  are 
shown  in  panel  (b),  one  without  added  noise,  and  one  with  noise  added.  The  error  bars  reflect  the 
amplitude  of  the  noise,  which  is  independent  Gaussian  random  noise  with  a  standard  deviation  a 
—  0.04.  To  our  knowledge,  there  is  nothing  that  prevents  correlated  measurement  errors  in  radio 
occultation  experiments  (see  below).  However,  the  results  here  are  easily  generalized  for  the  case 
of  nonuniform  and  correlated  errors. 

Panel  (c)  shows  the  results  of  model  inversion  using  the  linear  transformation  in  (5).  Results  for 
the  noisy  and  noise-free  cases  are  shown  with  solid  and  dashed  lines,  respectively.  The  noise-free 
inversion  is  indistinguishable  from  the  original  density  profile.  Error  bars  for  the  noisy  case  were 
computed  using  the  standard  error  propagation  formula: 

Cm  =  &Cd(G  t)4  (9) 


The  problem  with  the  inverse  method  adopted  so  far  is  error  propagation.  The  model  inversion 
in  panel  (c)  of  Fig.  3  exhibits  large  fluctuations  and  wide  confidence  intervals  despite  the  modest 
amount  of  noise  added  to  the  simulated  data.  The  problem  grows  with  decreasing  range  to  the  point 
of  introducing  artificial  layers  and  negative-going  excursions.  The  negative  excursions  and  layers 
that  might  easily  be  confused  with  intermediate  or  sporadic  E  layers  in  nature  are  unacceptable 
artifacts  and  render  direct  linear  transformation  suboptimal. 

The  cause  of  the  problem  is  illustrated  by  panel  (d)  of  Fig.  3.  The  histograms  represent  the 
singular  values  of  G,  ranked  top  to  bottom  from  largest  to  smallest.  The  condition  number  of 
the  inverse  transform  can  be  calculated  from  the  ratio  of  the  largest  to  smallest  of  the  singular 
values.  That  ratio  is  approximately  n  in  this  problem.  While  that  is  not  large  enough  to  give  rise 
to  numerical  instability,  it  has  serious  consequences  for  error  amplification.  Also  shown  in  panel  (d) 
are  the  singular  vectors  with  the  two  largest  and  two  smallest  singular  values.  Large  (small)  singular 
values  correspond  to  unstructured  (highly  structured)  singular  vectors  with  most  of  their  content 
at  large  (small)  ranges.  Features  that  resemble  the  structured  singular  vectors  in  the  noisy  data 
will  be  amplified  much  more  than  features  that  resemble  the  unstructured  ones.  This  explains  the 
preponderance  of  noise  amplification  and  ragged  oscillations  in  the  model  inversion  at  low  ranges. 
The  behavior  is  an  inherent  feature  of  the  inverse  Abel  transform. 


Two  conventional  strategies  exist  for  mitigating  the  oscillatory  behavior  of  the  linear  inversion. 
The  first  is  to  truncate  the  singular  value  matrix,  regarding  as  zero  those  values  that  fall  below 
some  threshold.  The  corresponding,  problematic  singular  vectors  and  associated  noise  amplification 
behavior  are  thus  eliminated.  A  less  drastic  approach  is  to  filter  the  singular  values,  modifying  them 
artificially  so  that  they  do  not  assume  too  low  values.  Such  so-called  damped  SVD  methods  include 
the  method  of  Tikhonov  regularization.  Panel  (e)  in  Fig.  3  shows  a  regularized  version  of  the  model 
inversion  with  error  bars  recalculated  according  to  (9).  The  forward  and  inverse  models  in  (6)  were 
recalculated  after  dividing  the  singular  values  A*  by  the  filter  coefficients 


fi 


\a 

A*  +  Q2 


(10) 


12 


Amplitude  Regularized  model  MaxEnt  model 

Figure  3:  Example  occupation  profile  inversion  comparing  different  methods,  (a)  Simulated  “true” 
electron  density  profile  nc(r).  (b)  Phase  profile  <f>(s)  calculated  using  (4).  Simulated  phase  data 
with  (solid)  and  without  (dashed)  noise  are  shown,  (c)  Direct  linear  inversion  of  the  simulated  data 
using  (5).  Models  with  (solid)  and  without  (dashed)  noise  are  shown,  (d)  Singular  value  analysis 
of  the  linear  transformation  (see  text),  (e)  Linear  inversion  using  Tikhonov  regularization.  The 
dashed  line  is  the  true  model,  (f )  Results  obtained  using  maximum  entropy.  The  dashed  line  is  the 
true  model. 
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where  ot  is  the  regularization  parameter.  This  is  equivalent  to  finding  the  model  that  minimizes  a 
weighted  combination  of  the  L2  norms  of  the  noisy  data  prediction  error  Grh  —  d  and  the  model 
itself  m.  The  regularization  parameter  was  set  to  10  in  this  example  in  an  attempt  to  reduce 
artifacts  while  retaining  fine  structure.  While  the  results  are  superior  to  the  unregularized  model 
in  terms  of  having  smaller  uncertainties  and  fewer  sporadic  oscillations,  the  profile  in  (e)  is  clearly 
distorted,  the  peaks  being  rounded  and  the  valley  being  filled  in.  Note  that  while  a  non-negativity 
constraint  can  be  incorporated  into  regularization,  we  have  not  done  so  here. 

Another  candidate  for  the  prior  probability  of  a  model  is  its  entropy,  a  combinatorial  measure 
from  which  partition  functions  in  statistical  mechanics  are  constructed.  Entropy  is  also  a  measure 
of  information  content  [55].  A  model  high  in  entropy  has  a  high  probability  of  occurrence,  all 
other  things  being  equal,  and  represents  minimal  independent  information,  making  it  minimally 
committal  to  unmeasured  data.  A  preference  for  high-entropy  models  suppresses  artifacts,  since 
every  feature  in  the  model  must  be  supported  by  the  data. 

In  the  maximum  entropy  problem,  the  prior  probability  is  made  to  be  a  function  of  the  en¬ 
tropy  (usually  es).  This  step  can  be  somewhat  controversial,  as  entropy  lies  outside  the  narrow, 
frequentist  definition  of  probability,  having  more  to  do  with  the  prior  state  of  knowledge  and  the 
plausibility  of  the  model.  As  [17]  points  out,  frequencies  are  difficult  to  obtain  in  scientific  mea¬ 
surements,  particularly  exploratory  ones  like  those  of  interest  here.  However,  [33]  and  others  have 
argued  that  Bayes  theorem  applies  equally  well  to  theoretical  probabilities  based  on  information 
theory  as  to  frequencies.  Numerical  simulations  offer  a  means  of  testing  the  argument. 

We  have  performed  the  occultation  profile  inversion  using  the  algorithm  described  recently  by 
[29].  Results  are  shown  in  panel  (f)  of  Fig.  3.  The  inversion  appears  to  produce  an  unbiased  repro¬ 
duction  of  the  original  density  profile,  even  in  the  vicinity  of  sharp  peaks  and  steep  gradients.  The 
recovered  density  approaches  but  never  crosses  zero  at  the  layer  boundaries.  The  overall  behavior 
is  superior  to  the  Tikhonov  regularization  solution  for  any  choice  of  regularization  parameter  a. 
Evidently,  the  noisy  data  contained  sufficient  information  to  specify  the  presence  of  steep  peaks 
and  shallow  valleys  in  the  profile.  The  constraint  that  the  model  be  positive  definite  and  in  some 
sense  smooth  was  sufficient  to  suppress  spurious  layers  and  other  artifacts  without  any  artificial 
truncation  of  the  model  space.  Similar  results  are  obtained  for  different  profile  shapes,  noise  levels, 
and  model  sizes. 

The  purpose  of  this  study  was  to  illustrate  some  shortcomings  of  standard,  direct,  linear  inver¬ 
sion  of  radio  occultation  data  and  to  compare  and  contrast  two  common  alternative  approaches, 
Tikhonov  regularization  and  maximum  entropy.  Both  reduce  the  production  of  spurious  layers  and 
artifacts  at  low  altitudes  inherent  in  the  structure  of  the  inverse  Abel  transform.  The  former  is  a 
fixed  linear  filter  that  works  by  attenuating  singular  modes  which  tend  to  amplify  noise.  The  danger 
with  this  approach  is  that  it  necessarily  removes  fine  structure  and  distorts  the  model  profile.  The 
latter  is  an  adaptive  nonlinear  filter  which  suppresses  features  in  the  model  that  are  not  supported 
in  the  data.  Maximum  entropy  produces  positive  definite  models  which  retain  their  fine  structure 
and  exhibit  little  distortion. 

It  is  imperative  that  direct  linear  inversion  of  radio  occultation  profiles  not  be  part  of  any  oper¬ 
ational  use  of  data  from  the  COSMIC  or  other  satellite  mission.  Artifacts  in  the  standard  inversion 
scheme  take  the  form  of  convincing  layers  at  low  altitudes  which  will  certainly  be  misinterpreted 
as  sporadic  E,  intermediate,  or  bottomside  layers.  This  problem  can  be  avoided  by  incorporating 
regularization  or,  better  still,  by  using  a  Bayesian  inversion  strategy  based  on  an  entropy  metric. 
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